POLYNOMIAL HOMOTOPIES FOR DENSE, SPARSE 
AND DETERMINANTAL SYSTEMS 



Os 
Os 
OS 



3 

Hi 

(N 



JAN VERSCHELDE 

Abstract. Numerical homotopy continuation methods for three classes of polynomial sys- 
tems are presented. For a generic instance of the class, every path leads to a solution and the 
homotopy is optimal. The counting of the roots mirrors the resolution of a generic system 
that is used to start up the deformations. Software and applications are discussed. 
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1. Introduction 

Solving polynomial systems numerically means computing approximations to all isolated 
solutions. Homotopy continuation methods provide paths to approximate solutions. The 
idea is to break up the original system into simpler problems. To solve the original system, 
the solutions of the simpler systems are deformed into the solutions of the original problem. 

Date: July 11, 1999. 

Research at MSRI is supported in part by NSF grant DMS-9701755, benefited from a post-doctoral 
fellowship at MSRI and is also supported in part by NSF grant DMS-9804846 at MSU. 
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This paper presents optimal homotopies for three different classes of polynomial systems. 
Optimal means that for generic instances of the classes there are no diverging solution paths, 
whence the amount of computational work is linear in the number of solutions. In the next 
section we list the principal key words, definitions and main theorems for dense, sparse and 
determinantal polynomial systems. The proofs of these theorems follow from the correctness 
of the homotopies. 



Path-following methods are standard numerical techniques (|3|, [|, | 7T| , ||123| , |125|| ) to 



achieve global convergence when solving nonlinear systems. For polynomial systems we can 
reach all isolated solutions. In the third section we describe the paradigm of Cheater's homo- 
topy or coefficient-parameter polynomial continuation ( |f75| , [f75| ). This paradigm 

allows to construct homotopies for which singularities only occur at the end of the paths. To 
deal with components of solutions we use an embedding method that leads to generic points 
on each component. This method is essential to numerical algebraic geometry |UH . 



From ESI we cite: "Algebraic geometry studies the delicate balance between the geomet- 



rically plausible and the algebraically possible". By a choice of coordinates we set up an 
algebraic formulation for a geometric problem that is then solved by automatic computa- 
tions. While this approach is extremely powerful, we might get trapped into tedious wasted 
computations after loosing the original geometric meaning of the problem. In section four we 
stress the geometric intuition of homotopy methods. Compactifications and homogeneous 
coordinates provide us the tools to generate the numerically most favorable representations 
for the solutions to our problem. In section five we arrive at the heart of modern homotopy 
methods where we outline specific algorithms to implement the root countsQ The counting 
of the roots mirrors the resolution of a system in generic position that is used as starting 
point in the deformations. 

Polyhedral methods occupy the central part of current research, as they are responsible for 
a computational breakthrough in numerical general-purpose solvers for polynomial systems. 
Section six is devoted to numerical software with an emphasis on the structure of the package 
PHC, developed by the author during the past decade. Another novel and exciting research 
development concerns the numerical Schubert calculus, which is one of the major new features 
in the second public release of PHC. The author has gathered more than one hundred 
polynomial systems that arose in various application fields. This collection serves as a test 
suite for software and a gallery to demonstrate the importance of polynomial systems to 
mathematical modelling. In section seven we sample some interesting cases. 

The reference list contains a compilation of the most relevant technical contributions to 
polynomial homotopy continuation. Besides those we want to point at some other works 
in the literature that are of special interest. Some user-friendly introductions to algebraic 



geometry appeared in recent years: see PJ, P?|, pTfl, with computational aspects in JT5] 



and []16l. As Newton polytopes have become extremely important, we recommend [ 130 1 and 



the handbook chapters See also [TIM] for the interplay between the combinatorics of 



polytopes and the (real) roots of polynomials. A recent survey that also covers polyhedral 
homotopies along with other polynomial continuation methods appeared in |>1] . 



lr The term "root count" was coined by Canny and Rojas [ |l3| while introducing mixed volumes to compu- 
tational algebraic geometry. 



polynomial homotopies for dense, sparse and determinantal systems 3 
2. Three Classes of Polynomial Systems 



The classification in Table [I] is inspired by f4"4" |. The dense class is closest to the common 
algebraic description, whereas the determinantal systems arise in enumerative geometry. 



system 


model 


theory 


space 


dense 


highest degrees 


Bezout 


pn 


projective 


sparse 


Newton polytopes 


Bernshtem 


(C*) n 


toric 


determinantal 


localization posets 


Schubert 


G mr 


Grassmannian 



Table 1. Key words of the three classes of polynomial systems. 



For the vector of unknowns x = (xi, x 2 , . . . , x n ) and exponents a = (ai, 02, . . . , a n ) G N n , 
denote x a = x^x^ 2 ■ ■ ■ x°£ . A polynomial system P(x) = is given by P — {p\,P2, ■ ■ ■ ,Pn), 
a tuple of polynomials p^ 6 C[x], i = 1, 2, . . . , n. 

The complexity of a dense polynomial p is measured by its degree d: 

p( x ) = Yl c * xa > d = deg(p), (1) 

0<ai+a 2 H ha n <d 

where at least one monomial of degree d should have a nonzero coefficient. The total degree 
D of a dense system P is D = Yli=i deg(pi). 

Theorem 2.1. (Bezout [[LSII j The system P(x) = has no more than D isolated solutions, 
counted with multiplicities. 

Consider for example 

P(x h x 2 ) = I 3 Xl + XlX \ I \ = ° n with total degree D = 4 x 4 = 16. (2) 

Although D = 16, this system has only eight solutions because of its sparse structure. 

The support A of a sparse polynomial p collects all exponents of those monomials whose 
coefficients are nonzero. Since we allow negative exponents (a e Z n ), we restrict x 6 (C*) n , 
C* = C\{0}. 

p(x) = ^ c a x a , Va G A : c a ^ 0, A C Z n , #A < cx). (3) 

The Newton polytope Q of p is the convex hull of the support A of p. We model the structure 
of a sparse system P by a tuple of Newton polytopes Q = (Qi,Q2, ■ ■ ■ ,Qn), spanned by 
A = (Ai,A 2 , . . . , An), the so-called supports of P. 

The volume of a positive linear combination of polytopes is a homogeneous polynomial in 
the multiplication factors. The coefficients are mixed volumes. For instance, for (Q 1 , Q 2 ), we 
write: 

2!vol 2 (A 1 Q 1 + X 2 Q 2 ) = V 2 {Q X , Q x )\\ + 2 • V 2 {Q X , Q 2 )X 1 X 2 + V 2 (Q 2 , Q 2 )\l (4) 

normalizing V 2 (Q, Q) = 2\vo\ 2 (Q). For the Newton polytopes of the system (g): 2!vol 2 (Ai(5i+ 
X 2 Q 2 ) = 4Af + 2 • 8A1A2 + 5AI. To interpret this we look at Figure [I] and see that multiplying 
Pi and P 2 respectively by Ai and A2 changes their areas respectively with A^ and X 2 . The 
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cells in the subdivision of Q\ + Q 2 whose area is scaled by AiA 2 contribute to the mixed 
volume. So, for the example in (0), the root count is eight. 




Qi Qi Qi + Qi 



Figure 1. Newton polytopes Qi, Q2, a mixed subdivision of Q\ + Q2 with volumes. 



Theorem 2.2. (Bernshtein JQ) A system P(x) = with Newton polytopes Q has no more 
than V n (Q) isolated solutions in (C*) n , counted with multiplicities. 



The mixed volume was nicknamed [JT3[] as the BKK bound to honor Bernshtein [[I]], Kush- 
nirenko [^TJ, and Khovanskii 



For the third class of polynomial systems we consider a matrix [C|X] where C G £,( m+r *) xm 
and X E ^(m+r)xr reS p ec tively collect the coefficients and indeterminates. Laplace expan- 
sion of the maximal minors of [C|X] in m-by-m and r-by-r minors yields a determinantal 
polynomial 

p(x)= £ sign(I, J)C[I]X[J], U = {l,2,...,m + r}, (5) 

ju j = u 
in J = 

where the summation runs over all distinct choices I of m elements of U. The partition 
{/, J} of U defines the permutation U 1— > (I, J) with sign (J, J) its sign. The symbols C[I] 
and X[I] respectively represent coefficient minors and minors of indeterminates. Note that 
for more general intersection conditions, the matrices [C|X] are not necessarily square. 

The vanishing of a polynomial as in (|5|) expresses the condition that the r-plane X meets 
a given m-plane nontrivially. The counting and finding of all figures that satisfy certain 
geometric conditions is the central theme of enumerative geometry. For example, consider 
the following. 



Theorem 2.3. (Schubert Let m,r > 2. In C m+r there are 

l!2!3!---(r-2)!(r-l)!-(mr)! 
m ' r ~ ml (m+1)! (m+2)! ■ • ■ (m+r-1)! ^ ' 

r-planes that nontrivially meet mr given m-planes in general position. 



This root count d m ^ r is sharp compared to other root counts, see [^] and ||114|| for examples. 



We can picture the simplest case, using the fact that 2-planes in C 4 represent lines in P 3 . 
In Figure 2 the positive real projective 3-space corresponds to the interior of the tetrahedron. 



POLYNOMIAL HOMOTOPIES FOR DENSE, SPARSE AND DETERMINANTAL SYSTEMS 



5 



Figure 2: m = 2 = r. 



In Figure 2, we see two thick lines meeting four given 
skew fine lines in a point. When not all input planes have 
the same dimension, but when the number of solutions is 
still finite, Pieri's formula provides a root count [SB, 

0- 



In [44] the problem is solved in chains of nested sub- 
spaces, using a cellation of the Grassmannian G mr of 
r-planes in C m+r . A localization poset models |56j the 
specialization of the solution r-plane when the input is 
specialized. 



Algorithmic proofs for the above theorems consist in two steps. First we show how to 
construct a generic start system that has exactly as many regular solutions as the root 
count. Then we set up a homotopy for which all isolated solutions of any particular target 
system lie at the end of some solution path originating at some solution of the constructed 
start system. 



3. The Principles of Polynomial Homotopy Continuation Methods 

Homotopy continuation methods operate in two stages. Firstly, homotopy methods exploit 
the structure of P to find a root count and to construct a start system p(°\x) = that has 
exactly as many regular solutions as the root count. This start system is embedded in the 
homotopy 

(x, t) = 7 (1 - t)F {0) (x) + tP(x) = 0, te [0, 1], (7) 

with 7 G C a random number. In the second stage, as t moves from to 1, numerical 
continuation methods trace the paths that originate at the solutions of the start system 
towards the solutions of the target system. 

The good properties we expect from a homotopy H(x.,t) = are (borrowed from [p4|): 

1. (triviality) The solutions for t = are trivial to find. 

2. (smoothness) No singularities along the solution paths occur. 

3. (accessibility) All isolated solutions can be reached. 



Continuation or path- following methods are standard numerical techniques (@, f|, ||, [ffl 



123| , 125|| ) to trace the solution paths defined by the homotopy using predictor- corrector 
methods. The smoothness property of complex polynomial homotopies implies that paths 
never turn back, so that during correction the parameter t stays fixed, which simplifies the 
set up of path trackers. A pseudo-code description of a path tracker is in Algorithm 13] 



The predictor delivers at each step of the method a new value for the continuation param- 
eter and predicts an approximate solution of the corresponding new system in the homotopy. 
Figure |3| shows two common predictor schemes. The predicted approximate solution is ad- 
justed by applying Newton's method as corrector. The third ingredient in path-following 
methods is the adaptive step size control. The step length is determined to enforce quadratic 
convergence in the corrector to avoid path crossing. 
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Algorithm 3.1. Following one solution path by an increment-and-fix predictor-corrector 
method with an adaptive step size control strategy. 



Input: #(x, t), x* G C": #(x*, 0) = 0, 

e > 0, maxJt, maxsteps. 
Output: x*, success if ||i?(x*, 1)|| < e. 



homotopy and start solution 
accuracy and upper bounds 
approximate solution if success 



t:=Q; k := 0; 
h := max _step size; 
old J := t; olds* := x* 
previous jk* := x*; 
stop := false; 

while t < 1 and not stop loop 
t := min(l, t + h); 
x* := x* + h(x* — previous s*); 
Newton(/7(x, t), x*, e, max success); 
if success 

then a := mm(Expand(h),maxstep-size); 

previous jk* := olds*; 

oldJb := t; olds* := x*; 
else h := Shrink(h); 

t := o/dJ; x* := olds*; 

end if; 

fc := k+ I; 

stop := (h < minstepsize) or (A; > max steps); 
end loop; 

success := (||fT(x*, 1)|| < e). 



initialization 
step length 
back up values for t and x* 
previous approximate solution 
combines stopping criteria 

secant predictor for t 
secant predictor for x* 
correct with Newton's method 
step size control 
enlarge step length 
go further along path 
new back up values 
reduce step length 
step back and try again 

augment counter 
stopping criteria 

report success or failure 



Following all paths can be done sequentially, one path at a time, or in parallel, with for 
each solution path the same sequence of values of the continuation parameter. The sequential 
path-following method has the advantage that the low overhead of communication 0] makes 
it very suitable to run on multi-processor environments. Note that the memory requirements 
are optimal. 




Figure 3. The secant and tangent predictor with step length h. 
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To solve repeatedly a polynomial system with the same coefficient structure P(c, x) = 0, 
the homotopy ([?]) is applied with p(°) = P(c°,x) = a system with random coefficients c°. 
Solving P(c°,x) = is no longer trivial, so the name cheater's homotopy |57j is appropriate. 
A similar idea appeared in WE, [761. For coefficients given as functions of parameters, a 



refined version of cheater's homotopy in |59] avoids repeated evaluation of those functions 
during path following: 

H(x,t) = P((l - [t-t(l -t) 7 ]) c ° + (t-t(l -t) 7 )c,x) = 0, te[0,l],7GC. (8) 

In [^] it is proven that with (§) all isolated solutions of P(c, x) = can be reached and 
that singularities can only occur at the end of the paths. 

Typically, when using a cheater's homotopy, the computational effort spent towards the 
end of the paths often accounts for most of the work. The main numerical problem is then 
to distinguish irrelevant solutions at infinity from ill-conditioned but possibly meaningful 
solutions. End games f4"3f , [[77], [7J| [7j|, |95| provide several procedures to approximate 
the winding number of a path. Recently, Zeuthen's rule was applied in [|5(J to determine 
numerically the multiplicity of an isolated solution. Multi-precision facilities are useful for 
evaluation of residuals and root refinement for badly scaled solutions. 

In most applications, the polynomial systems have real coefficients and invite the use of real 
homotopies. In Jl I ] it was conjectured and proven in |6(| that generically, real homotopies 
contain no singular points other than a finite number of quadratic turning points. At those 
bifurcation points pairs of real solution paths become imaginary or conversely, complex 
conjugated solution paths join to yield two real solution paths. We refer to 0, [PSfl , p3 
and p0, |oTJ for a discussion of numerical techniques to deal with quadratic turning points. 



A remarkable application of real homotopies in the real world consists in the finding of the 
relevant parameters of a polynomial system to maximize the number of real roots, see 
for the 40 real solutions for the Stewart-Gough platform in mechanics. 



In [93] the use of homotopy continuation to deal with overdetermined and components of 
solutions is discussed. Geometrically one slices the components of solutions with as many 
random hyperplanes as the dimension of the components. The solutions to the original poly- 
nomial system augmented with these random linear equations for the hyperplanes are generic 
points of the components, constituting the main numerical data to study those components. 
In particular, the number of generic points one obtains by this slicing procedure equals the 
sum of the degrees over all top-dimensional components of solutions. 



To make the algorithms of |J3J more efficient, in [93|, the following embedding of the 
polynomial system P(x) = is proposed: 



Pi(x) + \z = 0, 



1,2, 



. n 



X/ C 3 X 3 + Z 







(9) 



where the Aj's and c,'s are random complex numbers. This embedding has the advantage 



over the algorithms in |93j that fewer solution paths diverge. Solutions to the system 
with z = lie on a component of solutions. By Bertini's theorem, all solutions with z ^ 
are regular. In [93 1, it is proven that those solutions can be used as start solutions to reach 
all isolated solutions of the original polynomial system P(x) = 0. 
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The embedding (^) is performed repeatedly in the routine 'Embed' in the algorithm (copied 
from IP11) below. 



Algorithm 3.2. Cascade of homotopies between embedded systems. 

Input: P, n. system with solutions in C n 

Output: (£i, Xi, Zifl = §. embeddings with solutions 

So '■= P; initialize embedding sequence 

for i from 1 up to n do slice and embed 

Si := Embed (Si- 1, Zj); Zi = new added variable 

end for; homotopy sequence starts 

Z n := Solve(£ n ); all roots are isolated, nonsingular, with z n ^ 

for i from n — 1 down to do countdown of dimensions 

u p , s I Si \ homotopy continuation 

H l+l :- tS l+l + (l-*J^ i+i J; t : 1 ^ to remove z l+1 

Xi := limits of solutions of -f/j+i 

as t — > with 2j = 0; on component 

Zi := H i+ i('x,Zi ^ 0,t = 0); not on component: these solutions 

are isolated and nonsingular 

end for. 

This embedding allows the efficient treatment of overdetermined systems and other non- 
proper intersections. By perturbing the added hyperplanes and extending the generic points 
by continuation, interpolation methods can lead to equations for the components. 

4. The Geometry of the Deformations 

Homotopy methods have an intuitive geometric interpretation. In this section we illustrate 
the geometry of the three types of moving into special position: product, toric, and Pieri 
deformations. These can be regarded as three applications of the principle of continuity or 
conservation of number in enumerative geometry. 

Product homotopies deform polynomial equations into products of linear equations. In 
Figure |] we see the line configuration at the start and the ellipse-parabola intersection in the 
end. Note that complex space is the natural space for deformations. The other two complex 
conjugated intersection points could not be displayed in Figure f|. 

The sparser a system, the easier it can be solved. In Figure |5] we illustrate the idea of 
making a system sparser by setting up a so-called polyhedral homotopy that reduces this 
particular system at t = to a linear system. The lower hull of the Newton polytope of 
this homotopy induces a triangulation, which is used to count the roots. In particular, 
every cell in the triangulation gives rise to a homotopy with as many paths to follow as 
the volume of the cell. The other root for the example in Figure |5] can be computed with a 
homotopy obtained from P by the substitution of variables x% <— Xito 1 and X2 <— x 2 t~ 1 . This 
transformation pushes the constant monomial up, so that at t = we have the nonconstant 
monomials in the start system to compute the other root. 

Figure |B| displays a special and a general configuration of four lines. The basis has been 
chosen such that two of the four input lines are spanned by standard basis vectors. To 
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Figure 4. Intersection of quadrics: a degenerate and a target configuration. 




Figure 5. Triangulation of the Newton polytope of P with polyhedral ho- 
motopy P. 



compute all lines that meet four given lines, one of the four given lines is moved into special 
position so that it intersects two other given lines, see the left of Figure || The solution 
lines must then originate at those two intersection points and reach to the other opposite 
line while meeting the line left in general position. 

The constructions above are in a sense [0 "heuristic proofs". With the general position 
assumption we cheat a bit, avoiding the hard problem of assigning multiplicities. Making 



this so-called [127] "method of degeneration" rigorous was an important development in 



algebraic geometry. 

To deal with solution paths diverging to ill-conditioned roots or to infinity we need to com- 
pactify our space. Instead of polynomials in n variables we consider homogeneous forms with 
coordinates subject to equivalence relations. While mathematically all coordinate choices are 
equivalent, we select the numerically most favorable representations of the solutions. 

The usual projective transformation consists in the change of variables Xi := j*, for % = 
1,2,... ,n, which leads to the homogeneous system P(z) = 0. To have as many equations 
as unknowns, we add to this system a random hyperplane. Except for an algebraic set of 
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ei ei 



Figure 6. In P 3 two thick lines meet four given lines L±, L2, £3, and L4 in a 
point. At the left we see a special configuration and the general configuration 
is at the right. 



the coefficients of this added hyperplane, all solution paths are guaranteed to stay inside 
C n+1 when homotopy continuation is applied. We refer to [Q for numerical techniques that 
dynamically restrict the computations to n dimensions. 

For sparse polynomial systems, we introduce as in ||116|| a new variable for every facet of 



the Newton polytopes. The advantage of this more involved compactification is based on the 
observation that when paths diverge to infinity certain coefficients of the polynomial system 
become dominant. With toric homogenization the added variables that become zero identify 
the faces of the Newton polytopes for the parts of the system that become dominant. This 



compactification works in conjunction with polyhedral end games |43| which are summarized 
in Section I573L 



The natural way to compactify G mr is to consider a multi-projective homogenization ac- 
cording to rows or columns of the matrix representations for the planes. In addition, we 
have that the planes are equivalent upon a linear change of basis. Choosing orthonormal 
matrices to represent the input planes leads to drastic improvements in the conditioning of 



the solution paths, see [46| and [p.!4j| for experimental data. 



5. Root Counts and Start Systems 

The main principle is that counting roots corresponds to solving start systems. Algorithms 
to illustrate this principle will be shown for little examples for the three classes of polynomial 

systems. 

For dense polynomial systems, the computation of generalized permanents model the 
resolution of linear-product start systems. The algorithms to compute mixed volumes lead to 
polyhedral homotopies to solve sparse polynomial systems. The localization posets describe 
the structure of the cellation of the Grassmannian used to set up the Pieri deformations. 

5.1. Dense Polynomials modeled by Highest Degrees. A polynomial in one variable 
has as many complex solutions as its degree. A linear system has either infinitely many 
solutions or exactly one isolated solution in projective space. By this analogy |k| we see 
that Bezout's theorem generalizes these last two statements: a polynomial system has either 
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infinitely many solutions or exactly as many isolated solutions in complex projective space 
as the total degree. 

As the above presentation of Bezout 's theorem suggests, the simplest cases are univariate 
and linear systems, which are used as start systems. For the example system (H), a start 
system p(°) (x) = based on the total degree D is given by two univariate quartic equations 
x\ — c\ = = x\ — C2, where c\ and c 2 are randomly chosen complex numbers. Note that 
the computation of D = 4 x 4 models the structure of the solutions of as four solutions 
for X\ crossed with four solutions for x 2 . 



The earliest approaches of this homotopy appear in JTJJ, |JT|, [32], and were further 
developed in [[52], [|70[] , ||129|| . The book [nj contains a very good introduction to the practice 



of solving polynomial system by homotopy continuation. Regularity results can be found 
54] and ||131|| . While this homotopy algorithm has a sound theoretical basis, the total 



in 



degree is a too crude upper bound for the number of affine roots to be useful in most 
applications. 



Multi-homogeneous homotopies were introduced in JF^, |73| and applied to various problems 
in mechanism design, see e.g. ||118| , |119fl . Similar are the random product homotopies |55], |56| , 
applying intersection theory in |58|| , but less suitable for automatic procedures. For our 
running example ([|), we follow the approach of multi-homogenization and we group the 
unknowns according to the partition Z = {{^i}, {^2}}- The corresponding degree matrix 
Mz has in its (i, j)-th entry the degree of the i-th polynomial in the variables of the j-th set 
of Z. The 2-homogeneous Bezout bound Bz is the permanent of Mz- 



P(x 1 ,x 2 ) = 

{ x\ + XiX 2 + 1 
} x\x2 + X\x\ + 1 



M z = M {{xiUx2}} B z = per (M z ) 
[ 41 1 =4x2+3x1 

~ 3 2 =11 



(10) 



The computation of the permanent follows the expansion for the determinant, except for the 
permanence of signs, as it corresponds to adding up the roots when solving the corresponding 
linear-product start system: 



p(0)( 3 



an 



lite -ft*) 

i=l i=l 
3 2 











(11) 



In most applications the grouping of variables follows from their meaning, e.g.: for eigenvalue 
problems Ax = Ax, Z = {{A}, {x\, x 2 , ■ ■ ■ , x n }}. Efficient permanent evaluations in conjunc- 
tion with exhaustive searching algorithms for finding an optimal grouping were developed 
In case the number of independent roots equals the Bezout bound, interpolation 



111 



120 



methods 1105 are useful. 



Partitioned linear-product start systems were developed in ||109|| elaborating the idea that 
several different partitions can be used for the polynomials in the system. Motivated by 
symmetric applications | |107[ , general linear-product start system were proposed in |106| . 
These start systems are based on a supporting set structure S which provides a more refined 
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model of the degree structure of a polynomial system. 

TTT1 (12) 



S 



{x 1 }{x 1 ,x 2 }{x 1 }{x 1 } 
{x 1 }{x 1 ,x 2 }{x 1 }{x 2 } 



To compute the bound formally, one collects all admissible n-tuples of sets, picking one set 
out of every row in the set structure. 

B s = #{({xi}, {xi, x 2 }), ({aci}, {x 2 }), ({xi, x 2 }, {xi}), 

({x 1 ,x 2 }, {x 1 ,x 2 }), ({x 1 ,x 2 }, {xi}), ({xi,x 2 }, {x 2 }), (13) 
({xi}, {x 2 }), ({xi, x 2 }, {x 2 })({xi}, {x 2 })({a;i}, {x 2 })} 

Each admissible pair corresponds to a linear system that leads to a solution of a generic start 
system: 

p(o)/ x \ _ f (xi + cu)(xi + c i2 x 2 + ci 3 )(xi + ci 4 )(xi + ci 5 ) = ^ 
\ (xi + c 2 i)(xi + c 22 x 2 + c 23 )(xi + c 24 )(x 2 + c 25 ) = 1 

This start system has Bs = 10 solutions. In ||106||, the following theorems were proven. 



Theorem 5.1. Except for a choice of coefficients belonging to an algebraic set, there are 
exactly Bs regular solutions to a random linear-product system based on the set structure S. 

The proof of the theorem consists in collecting the determinants that express the degen- 
eracy conditions. These determinants are polynomials in the coefficients and vanish at an 
algebraic set. 

Theorem 5.2. All isolated solutions to -P(x) = lie at the end of some solution path 
defined by a convex-linear homotopy originating at a solution of a random linear-product 
start system, based on a supporting set structure for P. 

The idea of the proof is to embed the homotopy into an appropriate projective space and 
to consider the projection of the discriminant variety as an algebraic set for the continuation 
parameter. See |k| for an alternative proof. 



A general approach to exploit product structures was developed in [jSC| . For systems 



whose polynomials are sums of products one may arrive at a much tighter bound replacing 
the products by one simple product. An efficient homotopy to solve the nine-point problem 
in mechanical design was obtained in this way. 



The complexity of this homotopy based on the total degree is addressed in |]T0[ where 
a-theory is applied to give bounds on the number of steps that is needed to trace the 
solution paths. A major result is that one can decide in polynomial time whether an average 
polynomial system has a solution. A similar analysis of Newton's method in multi-projective 
space was recently done in |T7| . 



While the above complexity results apply to random systems, the problem of automatically 
extracting and exploiting the degree structure of a polynomial system is a much harder prob- 
lem. Finding an optimal multi-homogeneous grouping essentially requires the enumeration 
of all partitions [ |120j| . With supporting set structures one may obtain a high success rate, 



see |k| for a efficient heuristic algorithm. Recent software extensions for finding optimal 
partitioned linear-product start systems are in [|12S[|. 
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5.2. Mixed Subdivisions of Newton Polytopes to compute Mixed Volumes. For (Q), 
we collect the exponent vectors of the system P in the supports A: 



P(x u x 2 ) = 

j xf + X\X 2 + 1 = 
\ x\x 2 + X\x\ + 1 = 



A=(A 1 ,A 2 ) 
A 1 = {(0,0), (1,1), (4,0)} 
A 2 = {(0,0), (1,2), (3,1)} 



(15) 



The supports Ai and A 2 span the respective Newton polytopes Q\ and Q 2 . 



The Cayley trick |33|, Proposition 1.7, page 274] is a method to rewrite a certain resultant 
as a discriminant of one single polynomial with additional variables. The polyhedral version 



of this trick as in ||102| , Lemma 5.2] is due to Bernd Sturmfels. It provides a one-to-one 
correspondence between the cells in a mixed subdivision and a triangulation of the so-called 
Cayley polytope spanned by the points of Ai embedded in a (2n — l)-dimensional space. 
See |45| for another application besides mixed- volume computation. As in [45[], Figure 7 
gives a "one-picture proof" of this trick, displaying the Cayley polytope for the supports A 
in (|T5|). Note that this construction provides a definition for mixed subdivisions. 

The Cayley polytope is spanned by the 
points in Ax, where each point of Ai is ex- 
tended with n — 1 zero coordinates, and the 
points in Ai where each point in Ai is ex- 
tended with the respective z-th standard ba- 
sis vector, for i = 1, 2, . . . , n — 1. 

Omitting the added coordinates of this Cay- 
ley embedding, every cell in a triangulation 
of the Cayley polytope is identified with a 
cell in a mixed subdivision of the original tu- 
ple of polytopes. We can see this identifi- 
cation geometrically when slicing the Cayley 
polytope with a hyperplane that separates 
the embedded polytopes. As in Figure 7, the 
slice contains AiQi + X 2 Q 2 and the cells of a 
mixed subdivision are cut out by the cells in 
a triangulation of the Cayley polytope. 




AiQi + X2Q2 



Figure 7 : Cayley polytope of Q\ and Q 2 . 



The Cayley trick was implemented in [|112j| as an application of the dynamic lifting algo- 
rithm to construct regular triangulations. This method calculates the volume polynomial (f|) 
completely. When one is only interested in the mixed volume, the method is only efficient 
when the supports do not differ much from each other. 



To compute only the mixed volume, the lift-and-prune approach was presented in p4 



using a primal model to prune in the tree of edge-edge combinations. This approach operates 
in two stages. First the polytopes are lifted by adding one coordinate to every point in the 
supports. In the second stage, one computes the facets of the lower hull of the Minkowski 
sum that are spanned by sums of edges. These facets constitute the mixed cells in a mixed 
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subdivision. On the supports A in ([15D, we consider the lifted supports 
A=(A 1 ,A 2 ) 



{(0,0,1), (1,1,0), (4,0,0)} 
{(0,0,0), (1,2,0), (3, 1,1)} 

The lower hulls of the lifted polytopes are displayed in Figure ||. 



A x 
A 2 



(16) 




Q1 + Q2 



Figure 8. Lifted polytopes Qi, Q2, and a regular mixed subdivision of Qi + Qi 



The two cells that contribute to the mixed volume are identified by inner normals a and 
(3 that satisfy systems of linear equations and inequalities: 



a = (0,0,1) 
J 4«i = «i + a 2 < 
I «i + 2a 2 = < 



1 

3CKi + «2 + 1 



P= (2,-1,1) 

j fa + fo 

I ft + 2ft 



1 < 4ft 

< 3ft + ft + l (17) 



These systems express that the cells correspond to facets spanned by the sum of two edges 
on the lower hulls of Q\ and Q 2 respectively. The lift-and-prune method with a dual version 
of the linear inequality constrains as in (|17|) was elaborated in [ 112|| , exploiting the fact 



that several polynomials can share the same Newton polytope (see jJ0|) and with dimension 
reductions. 

The geometric dual construction to Figure | is displayed in Figure pi 





' mm 



Q1+Q2 



Figure 9. On the left we see the projection of a regular mixed subdivision of 
Qi + Q2- On the right, we have the dual construction with complexes J\fy(Qi) 
and MKQ2) collecting the cones of all vectors normal to the edges on the 
lower hulls of Q\ and Q 2 respectively. The intersection of the cones contain 
the normals to the mixed cells. 
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As in ||40|| , we assume that there are r different Newton polytopes. Given a tuple of 
lifted point sets A = (Ai,A 2 , ... , A r ), any lifted cell C v of a regular subdivision can be 
characterized by its inner normal as 

C v = (MiA4-AA.). (is) 

Since conv(CV) = conv(X^ = i d v Ai) is a facet of the lower hull, the inner product (., v) attains 
its minimum over Aj at 9 v Aj, i.e., 

Va,be<9 v Ai: (a,v) = (b,v), i = l,2, ...,r, (19) 



Va G Aj \ a v A, Vb G 9 V A; : (a, v) > (b, v), 



1,2, 



r. 



(20) 



Algorithm [5J] (presented in ||112|| ) gives a way to compute all mixed cells by searching for 
feasible solutions to the constraints ( |l9|) and (^0|). The algorithm generates a tree of all 
possible combinations of fcj-faces, with feasibility tests to prune branches that do not lead to 
mixed cells. The order of enumeration is organized so that mixed cells which share some faces 
also share a part of the factorization work to be done to solve the system defined by (|19|). 

Algorithm 5.1. Pruning algorithm with shared factorizations subject to inequality con- 
straints: 

Input: (Ai_,A 2 , ■ ■ ■ ,A r ), 

k = (ki, &2, ■ ■ ■ j k r ), n = y~]^—i ki, 

(J-'l, J-2, ■ ■ ■ , Tr). 

Output: e w = {Ces LJ \ V n {C, k) > }. 



lifted point sets 
Ai appears ki times 
ki-faces of lower hull of conv(Aj) 
collection of lifted mixed cells 



At level i, 1 < i < r: 
DATA and INVARIANT CONDITIONS: 

i-l 

(Mi, «): Miv = ^ v n+1 = 0,K = J2kj 

i=i 

(M 2 , «): M 2 v > ^ -v n+1 > 
dim(M 2 ) = n — k 
ALGORITHM: 

for each Ci G Ti loop 
Triangulate(M 1 , k, Q); 
if Miv = ^ w n+ i = 
then Eliminate(Mi, M 2 , re, C*j, Aj); 
if M 2 v > ^ > 

then proceed to next level i + 1; 
end if; 

end if; 
end for. 
At level z = r: 
Compute v: MiV = 0; 
Merge the new cell C v with the list & w . 



equalities (|19D 
upper triangular up to row k 

inequalities ( |20| ) 
sti/i feasible and reduced 

enumerate over all ki-faces 
ensure invariant conditions 

test for feasibility w.r.t. ( |T9|) 
eliminate unknowns 

test for feasibility w.r.t. (|20|) 
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Note that (|20|) has to be weakened to > type inequalities in order to be able to compute 
also subdivisions that are not fine. This also explains the merge operation at the end. 
The feasibility tests in the algorithm allow an efficient computation of the mixed cells. 
The conditions (|l~9|) and ( [20"t) are verified incrementally. After choosing a &«-face Ci = 
{c i, c u , . . . , c fc -j} of Ai, linear programming is used to check whether (Ci, ... ,Ci) can lead 
to a mixed cell in the induced subdivision. 

We end this section with complexity results. The complexity of computing mixed volumes 



is proven jn|] to be #P-hard. This complexity class is typical for all enumerative prob- 
lems, since, unlike the class NP, there exists no algorithm that runs in polynomial time for 
arbitrary dimensions to verify that a guessed answer is correct. Although the current algo- 
rithmical practice suggests that computing mixed volumes is harder than computing volumes 
of polytopes (which is also known as a #P-hard problem pOfl), this is not the case from a 
complexity point of view as shown in pl| . In [25|] it is shown that the mixed volume V n (Q) 



is bounded from below by n\vo\ n (Q M ) , being the polytope of minimum volume in Q. 

5.3. Sparse Polynomial Systems solved by Polyhedral Homotopies. The simplest 
system in the polytope model that still has isolated solutions in (C*) n has exactly two terms in 
every equation. Polyhedral homotopies solve systems with random complex coefficients 



starting from sparser subsystems. For (0), the homotopy with supports A as in (TTBI) is 



P(x l} X 2 ,t) = \ l % X K^T lX l t l^ C fl~n Wlth Cl , C ' t eC*. (21) 

The exponents of t are the values of the lifting u applied to the supports. 

To find the start systems, we look at Figure ||, at the subdivision that is induced by this 
lifting process. The start systems have Newton polytopes spanned by one edge of the first 
and one edge of the second polytope. Since the two cells that contribute to the mixed volume 
are characterized by their inner normals a and (3 satisfying ([17|) we denote the start systems 
respectively by P a and P 13 . To compute start solutions, unimodular transformations make 
the system triangular as follows. After dividing the equations so that the constant term is 
present, we apply the substitution X\ = y 2 , x 2 = y{ y\ on P a as follows: 



P a (x) 

The substitution in 
elaborated as 



>"■*-' ry 



- 1 + d[ 



o 
o 



p-(x = y ^) = | 
is apparent in the notation (as used in 



yi + c'{ 
yi 2 vl + 4 



■? — 

1 2 

/">"■ -*■ rf'-' 



{yblf 



30-l-(-l) 

2/i 

l-0+2-(-l) 

yi 



yi- 1 - 1 - 3 



The exponents are calculated by the factorization VU = L: 



' 3 -1 " 




1" 




1 " 


1 2 




-1 3 




-2 7 



r v 



(y 



o 



U\V 



y\ ■ y% 

-2 7 



(22) 



rVU 



(23) 



(24) 
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Since det(U) = 1, the matrix U is called unimodular. 



The polyhedral homotopy (|21f) directly extends the solutions of P a to the target system. 
To obtain a homotopy starting at the solutions of P@, we substitute in (|2~1~D x\ <— Xit^ 1 , 
X2 X2t and clear out the lowest powers of t. This construction appeared in and 
provides an algorithmic proof of the following theorem. 

Theorem 5.3. (Bernshtein 0, Theorem A}) For a general choice of coefficients for P, the 
system .P(x) = has exactly as many regular solutions as its mixed volume V n (Q). 



The original algorithm Bernshtein used in his proof was implemented in ||110|| 



For the numerical stability of polyhedral continuation, it is important to have subdivisions 
induced by low lifting values, since those influence the power of the continuation parameter. 
In ||112|] explicit lower bounds on integer lifting values were derived, but unfortunately the 



dynamic lifting algorithm does not generalize that well [65| if one is only interested in the 



mixed cells of a mixed subdivision. A balancing method was proposed in [[3(J to improve the 
stability of homotopies induced by random floating-point lifting values. 

Once all solutions to a polynomial system with randomly generated coefficients are com- 
puted, we use cheater's homotopy to solve any specific system with the same Newton poly- 
topes. One could say that polyhedral homotopies have removed the cheating part. The main 
advantage of polyhedral methods is that the mixed volume is a much sharper root count in 
most applications, leading to fewer paths to trace. They also allow more flexibility to exploit 
symmetry as demonstrated in [ |108|| . 

In case the system has fewer isolated solutions than the mixed volume, we consider the 
face systems. Define the face of a polynomial p with support A as follows: 

p(x) = ^2 c aX a has faces <9 v p(x) = 2J CaX - a for v 7^ 0. (25) 

a£A aedvA 

For v^O, the corresponding face system of P is d v P = (d v pi, d v p 2 , ■ ■ ■ , d v p n ). 

Theorem 5.4. (Bernshtein |7], Theorem B]) Suppose V n (Q) > 0. Then, P(x) = has fewer 
than V n (Q) isolated solutions if and only i/<9 v P(x) = has a solution in (C*) n , for v/0. 

As is the case for our running example (0), the Newton polytopes may be in generic 
position such that for any nonzero choice of the coefficients, the system has exactly as many 
isolated solutions as the mixed volume. In practical applications however, how can we decide 
whether paths are really going towards infinity? Relying on the actual computed values is 
arbitrarily, because 10 4 is as far from infinity as 10 8 , so we need algebraic structural data to 
certify the divergence. 

In the polyhedral end game (|3| solution paths are represented by power series expansions: 

: ?r£ + °<»» tM1 , SK0 . (2 6 ) 

The winding number m is lower than or equal to the multiplicity of the solution. For a solu- 
tion diverging to infinity or to a zero- component solution we observe that Vi 7^ 0. According 
to Theorem |5T^j, this solution vanishes at a face system d v P (same v with components V{ as 
in (PB|)), certifying the divergence. 
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To check whether a solution path really diverges is equivalent to the test on the value 
for V{. A first-order approximation of V{ can be computed by 

log |Xi(>i)| - log \Xi(s )\ 



Vi + O(s ) 



(27) 



log(si) - log(so) 

with < si < Sq. The above formula assumes the correct value for m. To compute m, 
solution paths are sampled geometrically with ratio h as s k = h k / m s . The errors on the 
estimates for v; are 



01) - (log \xi(s k +i)\ - log \xi(s k+2 )\) 



(28) 
(29) 



(\og\xi(s k ) \ - \og\xi(s k ^ 
= Cl h klm s {l + O{h k/m )). 

(k) 

An estimate for m is derived from two consecutive errors e\ . Extrapolation improves this 
estimate. So, by an inexpensive side calculation at the end of the paths, we obtain important 
structural algebraic information about the system. 

Recall that V n (Q) count the roots in (C*) n . Using Newton polytopes to count affine roots 
(i.e.: in C n instead of (C*) n ) was proposed in [^] with the notion of shadowed polytopes 
obtained by the substitution Xi <— Xi + q for arbitrary constants Cj. To arrive at sharper 
bounds, it suffices (see [^2| and |]86| ) to add a random constant to every equation. Sta- 
ble mixed volumes |42| provide a generically sharp affine root count. The constructions 
in f^(J and |29| avoid the use of recursive liftings to compute stable mixed volumes. Further 



developments and generalizations can be found in [3?|. 



5.4. 



Determinantal Polynomials arising in Enumerative Geometry. Homotopies for 

The algorithms in the numerical 



solving problems in enumerative geometry appeared in |44 
Schubert calculus originated from questions in real enumerative geometry |)6, [97]] and have 
their main application to the pole placement problem Ul2 |, [[S3, 84}], |j88| , |90|| in control theory. 



theorems, avoiding the 



The enumerative problems are formalized in some "finiteness" 
explicit but involved (as in (^)) formulas for the root counts. 

Theorem 5.5. The number of r-planes meeting mr general m-planes in C m+r is a finite 
constant. 



The first homotopy presented in [44] uses a Grobner basis for the ideal that defines G mr , as 



is derived in ]|101|| . By Grobner bases questions concerning any polynomial system are solved 
by relation to monomial equations. Every Grobner basis defines a flat deformation, which 
preserves the structure of the solution set p2| . Geometrically, this type of deformation is 
used to collapse the solution set in projective space to the coordinate hyperplanes, or in the 
opposite direction, to extend the solutions of the monomial equations to those of the original 
system. The flat deformations that are obtained in this way are similar to toric deformations 
in the sense that one moves from the solutions of a subsystem to the solutions of the whole 
system. 



The Grobner homotopies of j±4| work in the synthetic Pliicker embedding, and need to 
take the large set of defining equations of G mr into account. When expanding the minors 
into local coordinates, these equations are automatically satisfied, which leads to a much 
smaller polynomial system. Consequently, the second type of homotopies of the so- 
called SAGBI homotopies are more efficient. Instead of an ideal, we now have a subalgebra 
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and work with a SAGBI basis, i.e.: the Subalgebra Analogue to a Grobner Basis for an Ideal. 
The term order selects the monomials on the diagonal as the dominant ones. This implies 
that in the flat deformation 



sec 



103|| for a general description) only the diagonal monomials 



remain at t = 
For m = 2 



0. 



the equations of the SAGBI homotopy in determinantal form are 



Pi(x) = det 



r ji) 


Si) 






c n 






&12 


c 21 


li) 

c 22 


^21^ 


%22 


c 31 


Jt) 

c 32 


1 





c (i) 


c {i) 





1 


c 41 


c 42 





1,2,3,4, 



(30) 



where the coefficients c^j are random complex constants. In expanding the minors of (|30"D, 
the lowest power of t is divided out, minor per minor. The system at t = is solved 
by polyhedral continuation. The system at t — 1 serves as start system in the cheater's 
homotopy to solve any system with particular real values for the coefficients cjy . Figure [U] 
outlines the structure of the general solver. 



Binomial 

Systems 



polyhedral 
homotopy 



Generic 
Complex 
System 



flat 



deform. 



Generic 
Complex 
Problem 



cheater's 
homotopy 



Specific 

Real 
Problem 



Figure 10. The SAGBI homotopy is at the center of the concatenation. 



SAGBI homotopies have been implemented ||114|| to verify some large instances of input 
planes for which it was conjectured that all solution planes would be real. We refer to |89] 
and |98| for related work on these conjectures. In [[59[] an asymptotic choice of inputs is 
generated for which all solutions are proven to be real. 

The third type of homotopies presented in |44| are the so-called Pieri homotopies. Since 
they are closer to the intrinsic geometric nature, they are applicable to a broader class of 
problems. In particular, we obtain an effective proof for the following. 

Theorem 5.6. The number of r-planes meeting n general (m + 1 — ki)-planes in C m+r , with 
k\ + &2 + • • • + k n = mr, is a finite constant. 



Note that when all ki — 1, we arrive at Theorem |575|. For general fcj ^ 1, we are not aware 
of any explicit formulas for the number of roots. 

Figure |ll] shows a part of a cellular decomposition of G22 with the determinantal equations. 
We specialize the pattern X that represents a solution line by setting some of its coordinates 
to zero. This specialization determines a specialization of the input lines as follows: take 
those basis vectors not indexed by rows of X where zeroes have been introduced. The special 
line Sx for this example is as in (|31"D spanned by the first and third standard basis vector. 

Figure [11] pictures patterns of the moving 2-planes in the Pieri homotopy algorithm for 
the case (m,r) = (2,2), see Figure |6|. The bottom matrix is the general representation of a 
solution that intersects already the two input lines spanned by the standard basis vectors. 
At the leaves of the tree by linear algebra operations we can intersect with a third input line. 
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xu 


" 








x n 


" 


det 


Sx 










= 


det 


Sx 


£21 










^32 
£42 . 






£32 




det[L 3 |X] = det 



211 


" 


^21 








^32 





£42 . 



[1 4] 



[2 3] 



[2 4] 



Figure 11. Part of a cellular decomposition of the Grassmannian of all 2- 
planes. At the right we have the short-hand notation with brackets. The 
bracket [2 4] contains the row indices to the lowest nonzero entries in X. 



Moving down the poset, we deform form the left configuration in Figure ^| to the general 
problem. 

Denote by L\ and L2 the lines already met by X. At the leaves of the tree in Figure 
we intersect with the fourth line L 4 . The special position for the third line L 3 is represented 
by the matrix Sx, which intersects any X with coordinates as at the leaves of the tree. In 
the homotopy H(X, t) = we deform the line spanned by the columns of Sx to line L 3 , for 
t = to * = I. 



S 



x 



1 



1 





H{X,t) 



det(L 4 \X) = 
det{{l-t)S x - 



tL 3 \X) = 



*e[o,i]. (31) 



Every solution X(t) of H(X(t),t) = intersects already three lines: L%, L2 and L3. At the 
end, for t — 1, X also meets the line L3 in a point. 

The homotopy (|3l|) deforms two solution lines, starting at patterns which have their row 
indices for the lowest nonzero entries respectively as in [1 4] and in [2 3]. The correctness of 

justifies the formal root count using the localization 



this homotopy (proven in [4|J and 
poset. This combinatorial root count proceeds in two stages. First we build up the poset, 
from the bottom up, diminishing the entries in the brackets under the restriction that the 
same entry never occurs twice or more. Secondly, we descend from the top of the poset, 
collecting and adding up the counts at the nodes in the poset. More examples and variations 
are in |46l . 



To solve the general intersection problem of Theorem |5.6| , the special (m + 1 — A;j)-planes 
lie in the intersection of special m-planes. In the construction of the poset one has to follow 
additional rules as to ensure a solution that meets the intersection of special m-planes. We 
refer to 146|] for details. 



The third enumerative problem we can solve is formalized as follows. 

Theorem 5.7. The number of all maps of degree q that produce r-planes in C m+r meeting 
mr + q(m + r) general m-planes at specified interpolation points is a finite constant. 



POLYNOMIAL HOMOTOPIES FOR DENSE, SPARSE AND DETERMINANTAL SYSTEMS 



21 



In |83|, [84] and ||122|1 explicit formulas are given for this finite constant along with other 
combinatorial identities. Following a hint of Frank Sottile (see also ||100| ) and reverse engi- 
neering on the root counts in |53 , Pieri homotopies were developed in [[46] whose correctness 



yields a proof for Theorem 5.7. 



The analogue to Figure 11 for maps of degree one into G22 is displayed in Figure 12 



det 



Sx 



™0 

x 21 






x} 2 s 

x 32* 

x 42 i 



det 



Sx 



™0 

x ll 


■ 


™0 
x 21 


x 22 


™0 
x 31 


™0 
x 32 





™0 

x 42 J 



det[L n \X(s,t)] = det 



L, 



x n 

„,0 
x 21 




x 12 s 
x 22^ 
x 32^ 

x\ 2 t 



X 



31 



[2 5] 



[3 4] 



[3 5] 



Figure 12. Part of a cellular decomposition of the Grassmannian of maps of 
degree 1 that produce 2-planes in projective 3-space. The bracket notation at 
the right corresponds to a matrix representation of the coefficients of the 
map X(s, t). 



To solve the problem in Theorem |5?7] we need a special position for the interpolation 
points. By moving those to infinity, the dominant monomials in the maps allow to re-use the 
same special m-planes, whose entries should be considered modulo m + r. The homotopy to 
satisfy the n-th intersection condition is: 



H{X(s,t),s,t) 



det(Li\X(s u ti)) = i = l,2, ... ,n-l 
det((l - t)S x + tL n \X(s, t)) = 

\s-l)(l-t) + {s-Sn)t = te[Q,l] 

Note that the continuation parameter t moves the interpolation point from infinity, at 
(1,0), to the specific value (s,t) = (s n , 1). 



(32) 
[s,t) = 



See ||100|| for information on the selection of the input planes so that all maps are real. 

As an example of another problem in enumerative geometry we mention the 27 lines on a 
cubic surface in 3-space. According to [|8T], this is one of the gems hidden in the rag-bag of 
projective geometry. In | J92[j , the 27 lines are determined by breaking up the cubic surface 
into three planes in a continuous way such that each intermediate position is nonsingular. 
It is shown that this continuous variation is also valid in the real field. 

6. Numerical Software for Solving Polynomial Systems 

In computer algebra one wants to compute exactly as long as possible and to defer the 
approximate calculations to the very last end. Exactly the opposite way is taken in homotopy 
methods: here we use floating-point arithmetic and only increase the precision when needed. 

Next we mention programs with special features for polynomial systems. See [||, Chapter 
VIII] for a list of available software for path following. HOMPACK [ 73| , |124| | is a general 
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continuation package with a polynomial driver. It has been parallelized || [361, extended 



with an end game f9*5| , and upgraded ||126|| to Fortran 90. POLSYS_PLP [ |128|| provides 
linear-product root counts to be used in conjunction with HOMPACK90. The Fortran code 
for CONSOL is contained in |7I| , Appendix 6]. The C-program pss |66| applies homotopy 
continuation with verification by a-theory. Pelican [[!9|, f|l[] implements in C the polyhedral 
methods of [i0|. Efficient Fortran software for polyhedral continuation with facilities to 
compute all affine roots is used in P9"| . The computation of mixed volumes with the C 
program mvlp |23|, is a crucial step for sparse resultants [117|. A distributed version has 
been created in 



PHC is written in Ada and originated during the doctoral research of the author [113 



Executable versions were first released at the PoSSo open workshop on software fill 
public release of the sources is described in 
organized along four stages of the solver. 



The 



|115| |. The package is organized as a tool box, 
Figure O presents the flow of the solver. The 
package is menu-driven and file-oriented. A general-purpose black-box solver is available. 



1. Preconditioning 



4. Validation 



V 



o Coefficient Scaling 
o Reduction of degrees 



2. Root Counting 



o Bezout : degrees 

o Bernshtein : polytopes 

o Schubert : SAGBI/Pieri 



=> 



\ \ o Refining of the roots 
o Analysis of condition #s 



3. Homotopy Continuation 



o Fix continuation parameters 
o Choose Predictor-Corrector 



Figure 13. The four stages in the flow of the PHC solver. 



The new second release of PHC uses Ada 95 concepts in the construction of the math- 
ematical library. It is developed with the freely available gnu-ada compiler (currently at 
version 3. lip) on various platforms. To run the software no compilation is needed, as bina- 
ries are available for Unix Workstations: running SUN Solaris and SGI Irix, and Pentium 
PCs: running Linux and Solaris. The portability of PHC is ensured by the gnu-ada compiler. 

Another main feature of the second release are the homotopy methods for the Schubert 
calculus. Implementing those homotopies was a matter of plugging in the equations and 
calling the path trackers. The third release of the package should offer a more compre- 
hensive environment to construct homotopies, providing an easier access to the two main 
computational engines: mixed-volume computation and polynomial continuation. 



7. The Database of Applications 

The polynomial systems in scientific and engineering models are a continuing source of 
open problems. Systems that come from academic questions are often conjectures providing 
computational evidence in a developing theory. In various engineering disciplines polynomial 
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systems represent a modeling problem, e.g.: a mechanical device. The origin of a polyno- 
mial system matters when the original problem formulation does not admit well-conditioned 
solutions. As a general method to deal with badly scaled systems to compute equilibria of 
chemical reaction systems, coefficient and equation scaling was developed in [|69| , see also |7T], 
Chapter 5] and |gf . 



The collection of test systems is organized as a database and available via the author's 
web pages, A good test example reveals properties of the solution method and has a mean- 
ingful application. Besides the algebraic formulation it contains the fields: title (meaningful 
description), references (problem source), root counts (Bezout bounds and mixed volume), 
and solution list. 

Instead of producing a huge list with an overview, we pick some important case studies. 



katsura-n: (magnetism problem |47j) The number of solutions equals the total degree 
D = 2 n , so the homotopy based on D is optimal to solve this problem. Because the 
constant term is missing in all except one equation, the system is an interesting test 
problem for affine polyhedral methods. 

camerals: (computer vision ||28|| ) The system models the displacement of a camera be- 
tween two positions in a static environment |j23| . The multi-homogeneous homotopy is 
optimal for this problem, requiring 20 solution paths to trace instead of D = 64. 

gamentwo: (totally mixed Nash equilibria for n players with two strategies [p7 > |5"8]| ) This 
is another instance where mult i- homogeneous homotopies are optimal. The number of 
solutions grows like n\e~ l as n — > oo. The largest system that is currently completely 
solvable is for n = 8 requiring 14,833 paths to trace. Situations exist for which all 
solutions are meaningful. 

cassou: (real algebraic geometry) This system illustrates the success story of polyhedral 
homotopies: the total degree equals 1,344, best known Bezout bound is 312 (see [Q), 
whereas the mixed volume gives 24. Still eight paths are diverging to infinity and 



polyhedral end games JE]] are needed to separate those diverging paths from the from 
the other finite ill-conditioned roots. 

cyclic-n: (Fourier transforms ||, [|) For n = 7, polyhedral homotopies are optimal, with 
all 924 paths leading to finite solutions. For n > 8, the mixed volume overestimates 



the number of roots and there are components of solutions. In |94| the degrees of the 
components were computed for n = 8,9. There are 34940 cyclic 10-roots, generated by 
1747 solutions. 



pole28sys: (pole placement problem [T2j) This system illustrates the efficiency of SAGBI 
homotopies for verifying a conjecture in real algebraic geometry [p8]| . With the input 
planes chosen to osculate a rational normal curve, an instance with all 1,430 solutions 
real and isolated was solved in [ 1 14 1 . The problem is relevant to control theory Ipcfl . 

stewgou40: (mechanism design [[L8|]) Whether the Stewart-Gough parallel platform in 
robotics could have all its 40 solutions real was a notorious open problem, until recently, 
as it was solved by numerical continuation methods |]I8"|| . The problem formulation 
in |Ti| is highly deficient: the mixed volume equals 1,536 whereas only 40 solution 
paths will converge. 
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We emphasize that we have optimal homotopies for three classes of polynomial systems, 
but not for all possible structures. Although one can solve a modelling problem by a black- 
box polynomial-system solver, knowing the origin of the problem leads in most cases to more 
favorable algebraic formulations that help the resolution of a polynomial system. To produce 
really meaningful solutions one often has to be close to the source of the problems and be 
able to interact with the people who formulate the polynomial systems. 



In closing this section we list some notable usages of PHC Charles Wampler |121| used a 
preliminary version of PHC to count the roots of various systems in mechanical design. Root 
counts for linear subspace intersections in the Schubert calculus were computed by Frank 



Sottile, see |)8| for various tables. A third example comes from computer graphics. To show 



that the 12 lines tangent to four given spheres can all be real, Thorsten Theobald used PHC, 
choosing appropriate parameters in the algebraic formulation set up by Cassiano Durand. 

8. Closing Remarks and Open Problems 

The three classes presented in this paper are by no means exhaustive, but give an idea of 
what can be done with homotopies to solve polynomial systems. The root counts constitute 
the theoretical backbone for general-purpose black-box solving. Yet, the homotopy methods 
are flexible enough to exploit a particular geometrical situation, with guaranteed optimal 
complexity when applied to generic instances. 

From algebraic geometry formal procedures based on intersection theory count the number 
of solutions to classes of polynomial systems. Examples are the theorems of Bezout, Bern- 
shtein and Schubert. For these situations we construct a start system and have a homotopy 
to deform the solutions to this start system to the solutions to any specific problem. There 
are many other cases for which one knows how to count but not how to deform and solve 
efficiently. Research in homotopy methods is aimed at turning the formal root counts into 
effective numerical methods. As open problem we can ask for a meta-homotopy method to 
connect formal root counting methods to solving generic systems and deformation proce- 
dures. 

In most applications, only the real solutions are important. Once we know an optimal 
homotopy to solve the problem in the complex case, we would like to know whether all 
solutions can be real and how the real solutions are distributed. The reality question appears 
for instance in the theory of totally mixed Nash equilibria and in the pole placement problem. 
Finding well-conditioned instances of fully real problems can be done by homotopy methods. 
The finding of 40 real solutions to the Stewart-Gough platform [TS|] is perhaps the most 
striking example. The question is to find an efficient procedure to deform from the complex 
case to the fully real case. 
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ulated the author's research. The author is deeply indebted to all his co-authors. The 
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